#include "borwein.h"

// Calculate terms that compound Pi
void *calculateTerms(void *term){
	int *pointer = NULL;
	int t;
	mpf_t tmp;		// temporary for some arithmetic operations

	// Inicializing temporary variable
	mpf_init(tmp);

	pointer = term;
	t = *pointer;

	switch(t){
		case 0: // Sequence A
			while(counta < ITERATIONS){
				mpf_init(pi_struct.a[counta]);

				// Doing this: a[i] = sqrt(a[i-1])/2 + 1/(2*sqrt(a[i-1]));
				mpf_sqrt(pi_struct.a[counta], pi_struct.a[counta-1]);
				mpf_div_ui(pi_struct.a[counta], pi_struct.a[counta], 2);
				mpf_sqrt(tmp, pi_struct.a[counta-1]);
				mpf_mul_ui(tmp, tmp, 2);
				mpf_ui_div(tmp, 1, tmp);
				mpf_add(pi_struct.a[counta], pi_struct.a[counta], tmp);

				pthread_mutex_lock(&mutexa);
				counta++;
				pthread_mutex_unlock(&mutexa);
			}
			break;
		case 1: // Sequence B
			while(countb < ITERATIONS){
				if(countb <  counta || counta >= ITERATIONS){
					mpf_init(pi_struct.b[countb]);

					// b[i] = sqrt(a[i-1])*(1 + b[i-1])/(a[i-1] + b[i-1]);
					mpf_sqrt(pi_struct.b[countb], pi_struct.a[countb-1]);
					mpf_add_ui(tmp, pi_struct.b[countb-1], 1);
					mpf_mul(pi_struct.b[countb], pi_struct.b[countb], tmp);
					mpf_add(tmp, pi_struct.a[countb-1], pi_struct.b[countb-1]);
					mpf_div(pi_struct.b[countb], pi_struct.b[countb], tmp);

					pthread_mutex_lock(&mutexb);
					countb++;
					pthread_mutex_unlock(&mutexb);
				}
			}
			break;
		case 2: // Sequence P
			while(countp < ITERATIONS){
				if( countp < countb - 2 || countb >= ITERATIONS){
					mpf_init(pi_struct.p[countp]);

					// p[i] = p[i-1]*b[i]*(1 + a[i])/(1 + b[i]);
					mpf_mul(pi_struct.p[countp], pi_struct.p[countp-1], pi_struct.b[countp]);
					mpf_add_ui(tmp, pi_struct.a[countp], 1);
					mpf_mul(pi_struct.p[countp], pi_struct.p[countp], tmp);
					mpf_add_ui(tmp, pi_struct.b[countp], 1);
					mpf_div(pi_struct.p[countp], pi_struct.p[countp], tmp);
					
					countp++;
				}
			}
			break;
		default:
			break;
	}

	// Clearing temporary variable
	mpf_clear(tmp);

	pthread_exit(term);
}

// Algoritmo Borwein sequencial para calculo de pi
void sequentialBorwein(void){
	long i;
	mpf_t a0, b0, p0, a1, b1, p1;
	mpf_t tmp;		// temporary for some arithmetic operations


    // DUVIDA > QUAL DEVE SER A PRECISÃO??
	mpf_set_default_prec(32);

	// Inicializing variables
    mpf_init(a0);
    mpf_init(b0);
    mpf_init(p0);
    mpf_init(a1);
    mpf_init(b1);
    mpf_init(p1);
    mpf_init(tmp);

	mpf_sqrt_ui(a0, 2);				// a[0] = sqrt(2);
	mpf_set_ui(b0, 0);				// b[0] = 0.0;
	mpf_add_ui(p0, a0, 2);			// p[0] = a[0] + 2;

	for(i=0; i<ITERATIONS; i++){

		// a[i+1] = sqrt(a[i])/2 + 1/(2*sqrt(a[i]));
		mpf_sqrt(a1, a0);
		mpf_div_ui(a1, a1, 2);
		mpf_sqrt(tmp, a0);
		mpf_mul_ui(tmp, tmp, 2);
		mpf_ui_div(tmp, 1, tmp);
		mpf_add(a1, a1, tmp);

		// b[i+1] = sqrt(a[i])*(1 + b[i])/(a[i] + b[i]);
		mpf_sqrt(b1, a0);
		mpf_add_ui(tmp, b0, 1);
		mpf_mul(b1, b1, tmp);
		mpf_add(tmp, a0, b0);
		mpf_div(b1, b1, tmp);

		// p[i+1] = p[i]*b[i+1]*(1 + a[i+1])/(1 + b[i+1]);
		mpf_mul(p1, p0, b1);
		mpf_add_ui(tmp, a1, 1);
		mpf_mul(p1, p1, tmp);
		mpf_add_ui(tmp, b1, 1);
		mpf_div(p1, p1, tmp);

		mpf_set(a0, a1);	// a[i] <- a[i+1];
		mpf_set(b0, b1);	// b[i] <- b[i+1];
		mpf_set(p0, p1);	// p[i] <- p[i+1];
	}

	// Printing Pi value
	gmp_printf("%.6Ff\n", p1);

	// Clearing variables
    mpf_clear(a0);
    mpf_clear(b0);
    mpf_clear(p0);
    mpf_clear(a1);
    mpf_clear(b1);
    mpf_clear(p1);
    mpf_clear(tmp);
}

//Algoritmo Borwein paralelo para calculo de pi
void concurrentBorwein(void){
	void *status;
	long i;
	const int t[3] = {0, 1, 2};

	// < ACHO QUE ISSO NAO EH NECESSARIO>
	// Inicializing MP pointer variables
	//mpf_init(pi_struct.a);
	//mpf_init(pi_struct.b);
	//mpf_init(pi_struct.p);


	// Allocating spaces to 'a', 'b' and 'p' sequences
	pi_struct.a = (mpf_t*)malloc(ITERATIONS*sizeof(mpf_t));
	pi_struct.b = (mpf_t*)malloc(ITERATIONS*sizeof(mpf_t));
	pi_struct.p = (mpf_t*)malloc(ITERATIONS*sizeof(mpf_t));

	// Initializing all MP element variables
	mpf_init(pi_struct.a[0]);
	mpf_init(pi_struct.b[0]);
	mpf_init(pi_struct.p[0]);

	// Setting initial values
	mpf_sqrt_ui(pi_struct.a[0], 2);					// a[0] = sqrt(2);
	mpf_set_ui(pi_struct.b[0], 0);					// b[0] = 0.0;
	mpf_add_ui(pi_struct.p[0], pi_struct.a[0], 2);	// p[0] = a[0] + 2;

	counta = 1;
	countb = 1;
	countp = 1;

	// Initializing semaphores
	pthread_mutex_init(&mutexa, NULL);
	pthread_mutex_init(&mutexb, NULL);

	// Creating and executing threads
	pthread_create(&thread[0], NULL, calculateTerms, (void *)&t[0]);
	pthread_create(&thread[1], NULL, calculateTerms, (void *)&t[1]);
	pthread_create(&thread[2], NULL, calculateTerms, (void *)&t[2]);

	// Waiting for threads
	pthread_join(thread[0], &status);
	pthread_join(thread[1], &status);
	pthread_join(thread[2], &status);

	// Destroying semaphores
	pthread_mutex_destroy(&mutexa);
	pthread_mutex_destroy(&mutexb);

	// Printing Pi value
	gmp_printf("%.6Ff\n", pi_struct.p[ITERATIONS-1]);

	// Clearing all elements
	for(i=0; i<ITERATIONS; i++){
		mpf_clear(pi_struct.a[i]);
		mpf_clear(pi_struct.b[i]);
		mpf_clear(pi_struct.p[i]);
	}

	// Freeing memory alocated ??
    free(pi_struct.a);
    free(pi_struct.b);
    free(pi_struct.p);

	pthread_exit(NULL);
}

